#!/bin/python
"""
This script demonstrates the core idea of molecular dynamics simulations
by taking a one-dimentional harmonic ossicilator as an example.

It is simple, yet all ingredients of MD are incorporated.

To run this code, simply invoke:
   python md.py

(C) LT Kong konglt(at)sjtu.edu.cn
"""

# initialization
x = 0.         # equilibrium/initial position
v = 0.5        # initial velocity
k = 1.         # sping constant
f = -k*x       # force
m = 1.         # mass
a = f/m        # acceleration

dt = 0.1       # time step
hdt = dt * 0.5 # half of time step

# setup
nstep = 1000   # total MD simulation steps
nout  = 1      # frequency to output info
Ep    = 0.5*k*x*x # potential energy
Ek    = 0.5*m*v*v # kinetic energy
Et    = Ep + Ek   # total internal energy

# Open file to write the trajectory and thermal info
fp = open("md.log", "w")
fp.write("#istep x v Ep Ek Et\n")

# Trajectory in xyz format
dump = open("md.xyz", "w")

# Main MD loop
for istep in range(nstep):
  # velocity verlet
  v = v + a * hdt
  x = x + v * dt
  f = -k * x
  a = f / m
  v = v + a * hdt

  # thermo and output
  if (istep % nout ==0):
    Ep = 0.5*k*x*x
    Ek = 0.5*m*v*v
    Et = Ep + Ek
    fp.write(f"{istep} {x} {v} {Ep} {Ek} {Et}\n")
    dump.write(f"1\nTimestep: {istep}\n")
    dump.write(f"X {x} 0. 0.\n")

# cleanup
fp.close()
dump.close()
